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Abstract 

This paper describes the Direct Fourier Permuation Algorithm, an effi- 
cient method of computing Bit Reversal of natural indices [1, 2, 3, ... , 2 k ] in a 
vectorial manner (k iterations) and also proposes the Vectorial Digit Reversal 
Algorithm, a natural generalization of Direct Fourier Permutation Algorithm 
that is enabled to compute the r-digit reversal of natural indices [1, 2, 3, ... , r k ] 
where r is an arbitrary radix. Matlab functions implementing these two algo- 
rithms and various test and comparative results are presented in this paper to 
support the idea of inclusion of these two algorithms in the next Matlab Signal 
Processing Toolbox official distribution package as much faster alternatives to 
current Matlab functions bitrev 'order and digitrev order . 
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1 Introduction 



Recurrence relations in the Danielson-Lanczos Lemma allow for the immediate 
implementation of an explicitly recursive function for FFT computation. Each time 
it calls itself, a 2 k point FFT computation is reduced to two 2 fc_1 FFT compu- 
tations. This direct approach (divide et impera and explicit backward recursion) 
allows for the FFT computation algorithm to be expressed by decimation-in-time 
implementation (1). Subdividing the initial vector a; up to a set of pairs on which 
the FFT computation is very simple is sometimes called buterfly operation, a name 
suggested by the graphical representation of the computation [1]. 

function X = mkFFT(x) 
% x - lx2 k line vector 

% X - Discrete Fourier Transform of x calculated through 2-radix Decimation 
% in Time Fast Fourier Transform with Explicit Backward Recursion. 

N — max(size(x)); 
ifN == 1 
X = x- 

else (1) 

oddind = 1 : 2 : N; xodd = x(oddind); 

evenind = 2 : 2 : N; xeven = x(evenind); 

O = mkFFT(xodd); 

E = mkF FT (xeven); 

For k = 1 : 1 : N/2 

X(k) = 0(k) + (exp(-2 *pi*i*(k-l) /N)) * E(k); 
X(N/2 + k) = O(k) - (exp(-2 * pi * i * (k - 1) / N)) * E(k); 

end 
end; 

Implicitly, during this computation, a permutation of the argument vector x 
takes place whenever the function calls itself. Consequently, a composed permu- 
tation of the initial vector argument x will occur up to and within the innermost 
call. This permutation will be referred further in this paper as Fourier Permutation 
and is sometimes named the Buterfly Permutation [1], or Bit Reversal Permutation 
[2]- [9] due to the fact that the reversed bit representation of the permuted index of 
the initial 0-based index (see the following example) is implicitly sorted in ascending 
order. 



Example 1: 

permuted index bit reprezentation 



000 

4 100 

2 010 

6 110 

1 001 

5 101 

3 011 

7 111 



reversed reprezentation 0-based index 



000 





001 


1 


010 


2 


011 


3 


100 


4 


101 


5 


110 


6 


111 


7 
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The Fourier Permutation of the 1-based index [1, 2, 3, 4, 5, 6, 7, 8] is [1, 5, 3, 7, 2, 6, 4, 8]. 
The stages leading to this permutation through the explicit calls in (1) are the follow- 
ing: [1, 2, 3, 4, 5, 6, 7, 8] -> [1, 3, 5, 7, 2, 4, 6, 8] -> [1, 5, 3, 7, 2, 6, 4, 8] -> [1, 5, 3, 7, 2, 6, 4, 8], 
where a left to right reading shows the permutations operated from the first to the 
last (the innermost) call. 



2 Additive Constants Method 

In what follows we aim to compute the Fourier Permutation in an iterative 
manner, by other methods than the 'Bit Reversal' algorithms proposed in [1] - [9]. 

Our primary intention is to formulate an algorithm for computing the function 
Yjv = mkFPerm(N), where Y/v is the Fourier Permutation of the vector of indices 
[1,2, . . . , N] and N = 2 k , with acceptable efficiency. The formulation of such an 
algorithm requires that a recurrence relation of first order should be found between 
the consecutive components of the resulting vector Ym- 



Vp e 1,N- 1 : Y N (p+ 1) = Y N (p) +C N (p); Y N (l) = 1; (2) 
where Cjv(1 : N) is an additive constants vector to be determined. 

Example 2: 

C 8 (l:7) = [ (+4) (-2) (+4) (-5) (+4) (-2) (+4) ] 

Y 8 (l : 8) = [ 1 5 3 7 2 6 4 8] 

The Fourier Permutation of indices 1,8 can be computed taking into account 
that each component of the resulting vector is the sum of the preceding component 
and a constant that depends on N, and that the first component always has the 
value 1. 

The quantity to be found in all the odd rank positions of vector Cn (+4 in the 
example under discussion) will be further referred to in this paper as the trivial 
additive constant of the Fourier Permutation, all the other being called nontrivial 
constants. 

Definition 1. Let N = 2 k ,k G N*. We shall call the aditive constants of the 
Fourier Permutation of indices [1, 2, 3, 2 fc ] the (2 k — 1) components of vector Cn 
thus constructed: 

i. C 2 i = 1; 

ii. Vfc G N*, k > 2 : C 2 k = [2C 2fc -i,-2 fe + 3,2C 2 *-i]; 

The basic properties of the additive constants corresponding to the Fourier Per- 
mutation of indices [1, 2, 3, 2 fe ] are given in the following proposition: 

PROPOSITION 1. If k > 2 and N = 2 fe , the vector of the additive constants 
of the Fourier Permutation (Cn ) has the following properties: 

i. All odd-rank components store the value of the trivial additive constant: 
Vp G l.JV/2 : C N (2p - 1) = 2 k ~ 1 ; 

ii. All non-trivial additive constants are negative: 
Vp e l.JV/2 : C N (2p) < 0; 
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iii. The lowest non-trivial additive constant splits the Cn vector into two equal 
vectors: 

C N (1 : 2 fe - 1 - 1) = C N (2 k - 1 + 1 : 2 k - 1); 

iv. The value of the minor nontrivial constant: 
C N (2 k ~ 1 ) = -2 fe + 3; 

v. All non-trivial additive constants, except the minor nontrivial constant, are 
even: 

Vp G 1, (N/2 - 1) : \C N (2p)\mod2 = 0; 

vi. The components that are symetrically placed in relation to the minor non- 
trivial constant are equal: 

Vp G 1, (N/2 - 1) : C N (2 k -p) = C N (p); 

vii. Forward recursion: at step (k + 1) each of the vectors formed with the compo- 
nents on the left and on the right of the minor non-trivial additive constant, 
respectively, is double the vector of constants computed at step k: 

if Aq = 2 k and N 2 = 2 fe+1 then: C Na (1 : Aq - 1) = 2C Nl = C Na (Aq + 1 : N 2 ); 

viii. Forward recursion between minor nontrivial constants: 

if Aq = 2 k and N 2 = 2 k+1 then: CW 2 (Aq) = 2C Nl (2 fe " 1 ); 

Consequence: If N — 2 k and k > 2 , then the number of unique non-trivial 
additive constants of the Fourier Permutation is (k — 1), and the value of the 
trivial additive constant is 2 fc ~ 1 . Consequently the computation of the Fourier 
Permutation of indices [1, 2, 3, 2 k ] is reduced to: 

• the computation of these k additive constants and their distribution in a 
template vector of length 2 k — 1. 

• the computation of each component of the resulting vector as the sum of the 
preceding component and the corresponding additive constant. 

3 Computing Fourier Permutation through 
Additive Constants Algorithm 

Due to property (PI. iii), in order to obtain the template vector of additive 
constants it suffices to determine its first 2 fe ~ 1 components, i.e. its first 2 k ~ 2 non- 
trivial additive constants (as all the others, i.e. all odd rank components store 
the trivial additive constant). According to the above considerations, the matlab 
function for generating the additive constants of the Fourier permutation and the 
function for computing the Fourier Permutation of indices [1, 2, 3, 2 k ] can be 
written as: 

function V = mkCAPF(N) 

% VI = intermediate vector of the first 2 fe ~ 1 non-trivial additive constants 
% A = 2 fc ,fc>2 
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% ct = trivial constant 
ct = N/2; % property (Pl.i) 
k = log2(N); 
VI = [-2,-5]; 

% the two non-trivial additive constants corresponding to case k=3 
for p = 4 : k 

VI = 2 * VI; % property (Pl.vii) 

c = max{size{V I)); 

VI = [VI, VI(1 : c - 1)]; % property (Pl.iii) 
VI = [VI, VI(c) - 3]; % properties (Pl.vii, Pl.viii) 
end; 

c = max(size(V "/)); 

VJ = [VJ, V/(l : c - 1)]; % property (Pl.iii) 

V = zeros(l, N - 1) + ct; % property (Pl.i) 
V(2 : 2 : AT - 1) = V/; 

function V = mkFPerm(N) 

V = zeros(l,N);V(l) = [1]; 
CiV = mkCAPF(N); 

tov p = 1 : N — I 

V(p+l) = V(p) + CN{p); 
end; 

The complexity of the computation of permuted indices depends on the multi- 
plication operations in the mkCAPF function and on the addition iterated in the 
mkFPerm function. Consequently the complexity of computing function mkFPerm 
is 0(Nlog2(N)) and may decrease if multiplications are excluded from the compu- 
tational mechanism and the number of iterations decreases, possibly to k. 

4 Direct Fourier Permutation Method 

PROPOSITION 2. If N = 2 k and k > 2, the vector of the additive constants of 
the Fourier Permutation, Cn — mkCAPF(N), has the following properties: 

i. The sum of all the additive constants of the Fourier Permutation is: 

N-l 

Y, C N (i) = 2 k -l; 
i=i 

ii. The sum of all non-trivial constants of the Fourier Permutation is: 

N/2 

^C JV (2*-l) = -(2 fe - 1 -l) 2 ; 

i=l 

iii. The sum of the first 2 fc ~ 1 constants of the Fourier Permutation is 1: 

N/2 

£cW(i) = l; 
i=i 
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iv. The sum of any 2 k 1 consecutive constants of the Fourier Permutation is 1: 

N/2+p-l 

VpGMV/2: C n(*) = 1; 

i—p 



Consequences: 

By construction, for any i G 1, 2 fc_1 the difference between the rank (i + 2 k ~ 1 ) 
component and the rank i component of vector Y is the sum of the 2 fc_1 consecutive 
additive constants starting with (and including) the rank i constant. 
Consequently: 

Y N (2 k - 1 + l:2 k )= Y N (l : 2 fe " 1 ) + 1; 

Moreover, using properties (Pl.vii) and (P2.iii,iv) applied to case N = 2 k ~ 1 , 
it follows that the sum of any 2 k ~ 2 consecutive constants selected from the first 
(2 fc_1 — 1) components of Cn is 2. However, by construction, for any i G l,2 fc ~ 2 , 
the difference between the rank (i + 2 k ~ 2 ) component and the rank i component of 
vector Yjv is the sum of the 2 fc ~ 2 consecutive additive constants starting with the 
rank i constant: 

Y N (2 k - 2 + 1 : 2 fe " 1 ) = Y N (1 : 2 k - 2 ) + 2; 

The sum of any 2 fc ~ 3 consecutive constants selected from the first (2 fc ~ 2 — 1) 
components of Cn is 2 2 , and therefore: 

Y N (2 k - 3 + 1 : 2 k - 2 ) = Y N (1 : 2 fe " 3 ) + 2 2 ; 

And the procedure can go on until the first component of vector Cn is reached 
by successive truncations like those above: 

Cjv(1) = 2 fe " 1 ; ^(2) = + 2 fe - x ; 

PROPOSITION 3. If N = 2 k and k > 2, the Fourier Permutation of the indices 
[1, . . . , N] , has the following property: 

Vp G 0,(fc- 1) : Y N (2 k - p - 1 + 1 : 2 fe - p ) = Y N (l : 2 k - p - 1 ) + 2 P 



5 Direct Fourier Permutation Algorithm 

(Implicit Bit Reversal) 

By applying the properties mentioned in Proposition 2 and their consequences, 
the generating function of the the Fourier Permutation of indices 1, 2 k can be rewrit- 
ten in k iterations, as follows: 

function V = dfp(b, N) 

% N — 2 k ,k > 0; 

% b is a natural number; 

% V is the Fourier Permutation of indices [b, b + 1, b + 2 k — 1] 
V=[b];p2 = N/2; 
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while p2 > 1 

V = [V, V + p2] ; % one vectorial addition and one memory reallocation of V 
p2 = p2/2; % update by one division 
end; 

Let V — dfp(0, N) . Taking into account that the decimal number obtained 
by reversing the fc-bit representation of the decimal number 2 fc ~ p is 2( fc ~ 1 )~( fe ~ p ) = 
2 P_1 , the function that returns the decimal values corresponding to the binary rep- 
resentations obtained by reversing the fc-bit representation of all of the components 
of V is the following: 

function W = mkBlORevKBit(N) 
% N = 2 k ,k > 0; 
W=[0];p2 = l; 
while p2 < N/2 

W = [W,W + p2]- 

p2=p2* 2; 
end; 

PROPOSITION 4. On each iteration within mkBlORevK Bit function, the 
intermediate result W is ascendently ordered. 
Consequence: 

The result W — mkBlORevK Bit(N) is ascendently ordered. As the length 
of W is N — 2 k and the values in W are distinct natural numbers corresponding 
to binary fc-bit representations, it follows that the maximal item in W is 2 k — 1. 
Consequently W = [0, 1, 2, 3, . . . , 2 fe - 1]. 

The above considerations allow for the formulation of the following theorem: 

THEOREM 1. (Correctness and Complexity of Direct Fourier Permutation 
Algorithm) If V = dfp(b,N), N = 2 k , then by implication the vector V meets the 
relation: 

V = b + R, (Matlab formalism), 

where R is a permutation of W, specifically the one that corresponds to the re- 
versed fc-bit representations of the components of W (i.e. R = dfp(0,N), or in 
other words, reversed fc-bit representation of R is sorted in ascending order, i.e. 
R = bitrevorder(W)). The arithmetical complexity of the computation of V is 
0(N + k-l) . 

Note: The above theorem proves that the Fourier Permutation is uniquely 
determined by its first component and by its most important property formulated 
as Proposition 3. This is because neither Proposition 3 nor Proposition 2 really 
depends on the first component of the index to be permuted ([0, 1,2,..., 2 fc_1 ] or 
[1, 2, . . . , 2 fc ] or [6, b + 1, b + 2, . . . , b + 2 fe ~ 1 ]). In other words, all the properties 
mentioned in this paper regarding both the Fourier Permutation and the set of 
additive constants of the Fourier Permutation, including the above theorem, arc 
independent of the particular choice of the initial index (languages like C uses 0- 
based indexing while Matlab uses 1-based indexing) . These are the reasons why the 
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formalism in the above theorem has been chosen to unify the descriptions of both 
cases mentioned in Example 1. 



6 Vectorial Digit Reversal Generalization 

As a natural generalization of the Direct Fourier Permutation Algorithm we 
propose the following algorithm that is enabled to compute the r-digit reversal of 
natural index [1,2,3,..., r k ] for arbitrary radices r in a vectorial manner: 

function V = vdigitrevorder(N, r) 
%N = r k 
crN = N/r; 

KV = crN*[0:r-l] + l; 

V = KV; 

crLenV = r; 

crStep = 1; 

while crLenV < N 

if crLenV == r crSte P 

crStep = crStep + 1; 

crN — crN jr; 

KV = V + crN; 

crLenKV = crLenV; 

crLenV — 2 * crLenV; 
else 

KV = KV + crN; 
crLenV — crLenV + crLenKV; 
end; 

V=[V,KV}; 
end; 

The idea of the above algorithm is to compute (at each step of the iteration) the 
current kernel vector KV and to concatenate its value at the end of the currently 
calculated partial result V using the following rule: 

While V is still a partial result (i.e. length of V < N): 

• if the current length of V is equal to a natural power of the radix r p then 

— update the current kernel vector increasing the number of its components 
(up to r p ) and increasing the components themselves: 

KV = V + r k -P- x 

— update the current partial result doubling the number of its components: 

V = [V,V + r k -P- 1 } 

• else (the current length of V is between r p and r p+1 ) 

— update the current kernel vector increasing all of its components: 
KV = KV + r k -P- 1 

— update the current partial result increasing the number of its components 
with r p : 

V = [V, KV] 
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Another form of the same algorithm, in which the similarity to the Direct 
Fourier Permutation Algorithm is obvious, is the following: 



Vectorial Digit Reversal: 

function V = vdro(N, r) 
V=[l];pr = N/r; 
while pr > 1 



function V = dfp(b, N) 
V=[b];p2 = N/2; 
while p2 > 1 



Direct Fourier Permutation: 



KV = V; 



for cont = 1 : r — 1 



V = [V, KV + cont * pr}; 
end; 

pr = pr/r; 
end; 



end; 



p2 



V 



[V,V + p2\; 



P 2/2- 



We prefer the above formulation of the Vectorial Digit Reversal Algorithm for 
two reasons: it is less redundant and it enables us to observe that these two func- 
tions vdro(N,r), dfp(b,N), have the same arithmetic complexity when r = 2 and 
N = 2 k . Also, during various tests performed by the author, the above formula- 
tion of the Direct Fourier Permutation Algorithm has been proved to be the fastest 
Matlab script for computing bit reversal permutation. On the other hand, both the 
present mathematical formulation of the bit reversal computation and the theoret- 
ical arithmetic complexity obtained in Theorem 1 suggest that the Direct Fourier 
Permutation Algorithm defines the minimal computational effort (minimal arith- 
metic complexity) for computing bit reversal permutation of an arbitrary-based 
index [b, b + 1, . . . , b + 2 k ~ 1 ] in natural arithmetic, unless a stronger property than 
Proposition 3 can be formulated. 

Despite the Matlab formalism that has been chosen in the description of both of 
the above algorithms, up to this point of the present work, there was no particular 
hypothesis (concerning some particular implementation or some particular compu- 
tational advantage that could have been gained by programming in one specific 
language, medium or platform) being pursued, and consequently, all the above re- 
sults are purely arithmetical. In the following section we will see how the particular 
result of Theorem 1 that concerns complexity can be refined by making the most 
of the speed of the vectorized Matlab calculation. 



In this section we aim to obtain experimentally determined time-complexity re- 
sults for the two Matlab script functions dfp and vdro, coded above. 

Test variables: the arithmetic complexity of the dfp function doesn't depend 
on the starting value (b) of the index [b, b + 1, . . . , b + 2 k ~ 1 ] and consequently the 
dfp function will be tested only against increasing values of the length of the index 
(N = 2 k ). According to the above considerations, almost identical time-complexity 
results are expected to be found for the dfp(b,2 k ) and vdro{2 k ,2) computations. 
The vdro function will be tested against the variable (k, r). 



7 Benchmark 



7.1 Benchmark considerations 



9 



Computing the medium execution times: the medium execution time for 
each test variable (k and (k, r) respectively) will be the average of execution times 
cumulated over a great number of repetitions. Each medium execution time thus 
calculated will obviously depend on the performance of the computer running the 
tests, but the nature of variation of the medium execution time against the test 
variables is a characteristic of the Matlab implementations of the algorithms in 
themselves and does not depend on the performance of the computer. On the other 
hand, it is all the more necessary to establish with accuracy the medium execution 
time, as the execution time intended to be estimated is shorter (cases of immedi- 
ate practical interest, where test variables are small). In these cases the number 
of repetitions will be the greatest and will decrease gradually as the value of the 
test variables increases, up to a value that guarantees that the resulting vector of 
medium execution times is a fair statistical reflection of reality. In establishing the 
number of repetitions corresponding to each individual test variable, as well as in 
establishing the minimal number of repetitions corresponding to high values of the 
test variables, we will consider to be a fair reflection of reality such a representa- 
tion of the medium execution times in which the first symptom of convergence is 
present, i.e. the curve of medium execution times is nearly a smooth, ascendent one. 

Limitations of the present digitrevorder Matlab function: there are three 
reasons for replacing the existing digitrevorder Matlab function within the Signal 
Processing Toolbox: 

i) the arbitrary radix r is limited to the integer range from 2 to 36. In the 
proposed implementation of the Vectorial Digit Reversal Algorithm (vdro function) 
there is no such limitation. 

ii) the computation of the vdro proposed function is much faster than that 
required to be done in the present digitrevorder function. 

iii) unexpected behavior of the digitrevorder function could be sometimes ob- 
tained (digitrevorder(\ : 81, 3) - for example). This is because of a data validation 
issue: the power k of the radix r is computed in the digitrevorder. m (present file) 
as being: 

radixpow = floor(loglO(N)/loglO(radixbase)) 

and this instruction sometimes fails to return the correct result of the calculus (see 
floor(logl0(27)/logl0(3)), floor(logl0(8l)/logl0(3)), floor(logl0(7 7 )/logl0(7)) 
for example) and consequently, that instruction must be replaced. 

7.2 Testing the proposed functions against the existing 
digitrevorder Matlab function 

This section assumes that the tested functions are called with the following 
syntaxes: dfp(l,2 k ), vdro(2 k , 2), digitrevorder(l : 2 fe ,2) where k takes several in- 
creasing natural values, unless other explicit syntax is present. We preferred to test 
digitrevorder instead of bitrevorder function because, due to its implementation, 
the latter is calling the former. 

As expected, the time-complexity lines (Figure 1) of the two proposed algo- 
rithms (dfp and vdro) are almost identical when the radix is r = 2. In this case 
the graphical representation in Figure 1 enables us to distinguish two variation 
regimes along the time-complexity lines of the both algorithms. Figure 2 reveals 
more clearly the point where the variation of the time-complexity line undergoes 
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Figure 1: Testing dfp, vdro, digitrev 'order against increasing k values 



Theoretical time-complexity lines 



-e- dfp (Direct Fourier Permutation) 

• Minimal theoretical time-complexity: Cmin*(N+k-1 ) 

. Maximal theoretical time-complexity: Cmax*(N+k-1 ) 

+ Maximal experimental time-complexity line (2 1 -2 12 ): Cmed*(4*k) 



8 10 12 14 

k=1 :20 (Array length: N=2 k ) 



20 



Figure 2: Time-complexity lines of the dfp computation 
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10" 



vdro versus digitrevorder: time-complexity lines 
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-O- digitrevorder(1 :5 ,5) 
O digitrevorder(1 :7 k ,7) 
-+- vdro(5 k ,5) 
vdro(7 k ,7) 




1.5 



2.5 3 3.5 4 4.5 
k=1 :6 (Array length: N=2 k ;) 



5.5 



Figure 3: Time-complexity lines: vdro versus digitrevorder (arbitrary radices) 



a significant change. During all the tests done by the author so far, the existence 
of this point has been proved to be an invariant of the algorithm. This is because 
the medium execution time calculated for each k depends on the following factors: 
k itself, the number of repetitions, the type of the cache and influences from other 
processes that are running on the same CPU. The greater the increase of k and in 
the number of repetitions, the greater the contribution of the two other factors will 
be to the calculated medium execution time and the computation itself will becomes 
less cache-optimal. In Figure 2, Cmin, Cmed and Cmax are the minimal, medium, 
maximal duration, respectively, of the computational cycle of the dfp algorithm, all 
of them being experimentally determined. 



Also, Figure 2 and the dfp function itself allow for the formulation of the fol- 
lowing remarks: 

i) due to its simplicity, the dfp function certainly has a lower complexity and 
a higher performance than the bitrevorder function within the Signal Processing 
Toolbox. 

ii) the dfp function allows us to compute the bit reversal order of indices 1, 2 k 
in just k iterations, each of these involving only three operations: one vectorial 
addition, one memory reallocation, and one division. 

iii) for k ranging between 2 and 12, the medium execution time needed to com- 
pute V = dfp(l, 2 k ) is increasing linearily with 4k (reflecting the very few operations 
within each iteration). On this range of k values, the dfp function takes all compu- 
tational advantage of its simplicity. Consequently, this is the range on which the dfp 
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function reaches its maximal efficiency and its minimal computational complexity 
0(4fc). 

iv) for k ranging between 13 and 20, due to the increasing size of variable V, 
the computational complexity of the dfp function suddenly turns into its own worst 
case scenario - predicted by its theoretical complexity 0(N + k, — I) - on which 
the medium execution time needed to compute V = dfp(l,2 k ) is increasing expo- 
nentially with k (i.e. linearily with N = 2 k ). But even on this range of k values, 
the multiplicative constant characterising the time-complexity of the computation 
is relatively small. Consequently, on this range of k values the dfp function is still 
operative, even though it reaches its minimal efficiency and its maximal computa- 
tional complexity 0(2 k + k — 1). 

Figure 3 reveals the great improvement in the performance of the proposed vdro 
function compared to the present Matlab function digitrevorder. Even for arbi- 
trary radices r the computation of the vdro function is faster and more reliable. 

7.3 Replicability of results 

The essential qualitative results illustrated by the tests above are the fol- 
lowing: 

i) higher performance of the dfp function compared to the Matlab bitrevorder 
function; 

ii) linear variation proportional with Ak of the medium execution time needed 
to compute V = dfp(l, 2 k ) for k ranging between 2 and 12. 

iii) exponential variation proportional with (2 k +k — l) of the medium execution 
time needed to compute V = dfp(l, 2 k ) for k ranging between 13 and 20. 

iv) higher performance of the vdro function compared to the Matlab digitrevorder 
function for arbitrary radices r; 

The replication of the qualitative results mentioned above does not depend on 
the PC hardware architecture used. Nonetheless, the replication of any quantitative 
results requires a hardware architecture which must be very similar to the one used 
for the purposes of this paper. All the tests mentioned above were run in Matlab 6.5 
R13, on the following configuration: CPU: Prescott, Intel Pentium 4E, 2800Mhz; 
Memory bus features: Type - Dual DDR SDRAM, Bus width - 64 bit, Real clock - 
200MHz (DDR), Effective clock - 400MHz, Bandwidth 6400MB/s; Memory module 
features: Size - 512 MB, Module type - Unbuffered, Speed - PC3200 (200MHz); 

As far as the variation of experimental medium execution times along the range 
of k is concerned, all the tests that have been run have high statistical relevance, 
being based on a number of repetitions that is high enough for statistic phenomena 
to become observable. All the above tests have been validated by the author by 
performing more tests on other hardware configurations, with similar results. 
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8 Appendix 



1. Discrete Fourier Transform: 
The Discrete Fourier Transform of the signal of finite length x(0 : N — 1) is vector 
X(0 : N — 1), with the following components: 

N-l 

Vfc G 0,(JV- 1) : AT(fc) = ^ x(n)W% n , W N = exp(-2ni/N); 

n=0 

2. Fast Fourier Transform, Danielson-Lanczos Lemma: 

Fast Fourier Transform is the algorithm that allows for the computation of the 
Discrete Fourier Transform with a complexity of order 0(Nlog2N). The FFT algo- 
rithm for computing the Discrete Fourier Transform of signals of length 2 k is based 
on the Danielson-Lanczos Lemma: 

Danielson-Lanczos Lemma: 
Let signal x(0 : N - 1), N = 2 k and /(0 : N/2 - 1), g(0 : N/2 - 1) be defined by: 

f(n) = x(2n), g(n) = x(2n + 1), n e 0, (N/2 - 1); 

Let X(0 : AT — 1), F(0 : AT/2 - 1), G(0 : A^/2 - 1) be the Fourier transforms of 
signals x, f, and g, respectively. 
Then X has the following components: 

X(k) = F(k) + W%G{k), k e 0, {N/ 2 - 1); 
AT(A72 + k) = F(k) - W^G(k), keO, (N/2 - 1); 
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